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Abstract. We present a multiscale approach to the modeling of polymer dynamics in the 
presence of a fluid solvent. The approach combines Langevin Molecular Dynamics (MD) techniques 
with a mesoscopic Lattice-Boltzmann (LB) method for the solvent dynamics. A unique feature of 
the present approach is that hydrodynamic interactions between the solute macromolecule and the 
aqueous solvent are handled explicitly, and yet in a computationally tractable way due to the dual 
particle-field nature of the LB solver. The suitability of the present LB-MD multiscale approach is 
demonstrated for the problem of polymer fast translocation through a nanopore. We also provide 
an interpretation of our results in the context of DNA translocation through a nanopore, a problem 
that has attracted much theoretical and experimental attention recently. 
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1. Introduction. Mathematical modeling and computer simulation of biolog- 
ical systems is in a stage of burgeoning growth. Advances in computer technology 
but also, perhaps more importantly, breakthroughs in simulational methods are help- 
ing to reduce the gap between quantitative models and actual biological behavior. 
The main challenge remains the wide and disparate range of spatio-temporal scales 
involved in the dynamical evolution of complex biological systems. In response to 
this challenge, various strategies have been developed recently which are in general 
referred to as "multiscale modeling" . Some representative examples include hybrid 
continuum- molecular dynamics algorithms pQ, heterogeneous multiscale methods [2], 
and the so-called equation- free approach [3]. These methods combine different lev- 
els of the statistical description of matter (for instance, continuum and atomistic) 
into a composite computational scheme, in which information is exchanged through 
appropriate hand-shaking regions between the scales. Vital to the success of this 
information exchange procedure is a careful design of proper hand-shaking interfaces. 

Kinetic theory lies naturally between the continuum and atomistic descriptions, 
and should therefore provide an ideal framework for the development of robust multi- 
scale methodologies. However, until recently, this approach has been hindered by the 
fact that the central equation of kinetic theory, that is, the Boltzmann equation, was 
perceived as an equally demanding approach as molecular dynamics from the com- 
putational point of view, and of very limited use for dense fluids due to the lack of 
many-body correlations. As a result, multiscale modeling of nanoflows has developed 
mostly in the direction of the continuum/molecular dynamics paradigm [1]. 

Over the last decade and a half, major developments in lattice kinetic theory [H[5] 
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are changing the scene. Minimal forms of the Boltzmann equation can be designed 
on the lattice, which quantitatively describe the behavior of fluid flows in a way that 
is often computationally more advantageous than the continuum approach based on 
the Navier-Stokes equations. Moreover, lattice kinetic theory has proven capable of 
dealing with complex flows, such as flows with phase transitions and strong hetero- 
geneities, for which continuum equations are exceedingly difficult to solve, if at all 
known (for a recent review see [6]). These advances have opened the road to devel- 
oping new mesoscopic multiscale solvers [7]. The present work provides a successful 
implementation of such an approach. We will focus on the coupling of a mesoscopic 
fluid solver, the lattice Boltzmann method, with simulations at the atomistic scale 
employing explicit molecular dynamics. A unique feature of our approach is the dual 
nature of the mesoscopic kinetic solver, which propagates coarse-grained informa- 
tion (the single-particle Boltzmann probability distribution), along straight particle 
trajectories. This dual field/particle nature greatly facilitates the coupling between 
the mesoscopic fluid and the atomistic levels, both on conceptual and computational 
grounds. 

The paper is organized as follows. In Section fj^l we present the basic elements 
of the multiscale methodology, namely the Lattice Boltzmann treatment of the fluid 
solvent, and its coupling to a Molecular Dynamics simulation of the solute biopolymer. 
In Section £J51 we present an application of this multiscale methodology to the problem 
of long polymer translocation through a nanopore; in particular, we analyze in detail 
the role of hydrodynamics in accelerating the translocation process. In Section SjUwe 
elaborate on the relevance of our results to the problem of DNA translocation, which 
has attracted much theoretical and experimental attention recently. We conclude in 
Section |j5]with general remarks and outlook for future extensions. 

2. Lattice-Boltzmann - Molecular-Dynamics multiscale methodology. 

We consider the generic problem of tracing the dynamic evolution of a polymer 
molecule interacting with a fluid solvent. This involves the simultaneous interaction 
of several physical mechanisms, often acting on widely separate temporal and spatial 
scales. Essentially, these interactions can be classified in three distinct categories as 
solute-solute, solvent-solvent and solvent-solute. The first category includes the con- 
servative many-body interactions among the single monomers in the polymer chain. 
Being atomistic in nature, these interactions usually set the shortest scale in the over- 
all multiscale process. They are typically handled by Molecular Dynamics techniques 
for constrained molecules. The second category, the solvent-solvent interactions, refer 
to the dynamics of the solvent molecules, which are usually dealt with by a continuum 
fluid-mechanics approach; in the present work these will be described by the meso- 
scopic Lattice Boltzmann equation. The second and third category have also been 
handled by simulating the solvent explicitly via molecular dynamics, implicit solvent 
particles via Brownian dynamics including hydrodynamic interactions, or solving the 
corresponding Fokker-Planck equation [8]. Finally, the solvent-solute dynamics will 
be treated by augmenting the molecular dynamics side with dissipative fluid-molecule 
interactions (Langevin picture) and including the corresponding reaction terms in the 
fluid-kinetic equations. 

2.1. Atomistic dynamics. We consider a polymer consisting of N monomer 
units (also referred to as beads). The polymer is advanced in time according to the 
following set of Molecular Dynamics-Langevin equations for the bead positions f p and 
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velocities v p : 



(2.2) = F p +F p J + F p + F p , p= 1,N 



where we distinguish four types of forces: 

(2.3) F p C =-Y^d Fp V(r p -r q ) 

9 

(2.4) F p =j(u p -v p ) 

(2.5) F p = m& 

(2.6) F p = -Xpdf p K p 

The first term represents the conservative bead-bead interactions through a potential 
which we will take to have the standard 6 — 12 Lennard- Jones form, 

(2.7) V LJ (r) = Ae[(a/r) 12 - (cr/r) 6 ] 

truncated at a distance of r=2 1 / 6 cr [9]. This was combined with a harmonic part to 
account for the energy cost of distorting the angular degrees of freedom, 



(2-8) V ang (</>) = ^~ 



2 



with <j) the relative angle between two consecutive bonds. Torsional motions will not 
be included in the present model, but can easily be incorporated if needed. 

We consider next the solute-solvent interactions. The second term on the right- 
hand-side of Eq. (|2.ip represents the mechanical friction between the single bead and 
the surrounding fluid, v p being the bead velocity and u p the fluid velocity evaluated 
at the bead position. In addition to mechanical drag, the polymer feels the effects 
of stochastic fluctuations of the fluid environment, through the random term, £ p , a 
Gaussian noise obeying the fluctuation-dissipation relations: 

< e P > = o 

< Ifo. t)fc q , t 1 ) > = ~f(k B T/m)V6(r p - r q )8{t - t 1 ) 

where V is the volume of the cell to which beads p and q belong. Finally, A p df p n p is 
the reaction force resulting from N — 1 holonomic constraints for molecules modelled 
with rigid covalent bonds: 

(2.9) k p = \r p+1 - r p \ 2 - rl = 

ro being the prescribed bond length, and {X p } is the set of N — 1 Lagrange multipli- 
ers conjugated to each constraint. The usage of constraints instead of flexible bond 
lengths makes it possible to eliminate unimportant high-frequency intra-molecular 
motion which would render the underlying LB propagation prone to numerical insta- 
bilities. In this way, the time-step of the Molecular Dynamics part can be increased by 
about one order of magnitude, as much as the overall efficiency of the LBMD method, 
as we shall discuss in section i j2.5l Finally, in order to avoid spurious dissipation, the 
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bead velocities are required to be strictly orthogonal to the relative displacements. 
Given the second order atomistic dynamics, the velocities must obey the independent 
constraints: 

(2.10) = (f p+1 - r p ) ■ (v p+1 -v p ) = Q 

The constraints (|2.9j) , (12.10p are enforced over positions and momenta separately via 
the SHAKE [ID] and the RATTLE algorithms [IT]- The implementation of these con- 
straints requires the iterative solution of the system of equations (|2.9|) - (|2.10[) . typically 
accomplished via standard Newton-Raphson techniques. 

2.2. Fluctuating Lattice Boltzmann method. The Lattice Boltzmann equa- 
tion is a minimal form of the Boltzmann kinetic equation in which all details of molecu- 
lar motion are removed except those that are strictly needed to recover hydrodynamic 
behavior at the macroscopic scale (mass-momentum and energy conservation). The 
result is an elegant equation for the discrete distribution function fi(x,t) describing 
the probability to find a particle at lattice site x at time t with speed v ' = C{. More 
specifically, since we are dealing with nanoscopic flows, in this work we shall consider 
the fluctuating Lattice Boltzmann equation which takes the following form: 

(2.11) fi(x + CtAt, t + At) = fi(x, t) - LoAt{U - f* q )(x, t) + F t At + S, At 

where fi (x, t) represents the probability of finding a fluid particle at spatial location 
x and time t with discrete speed ci. The particles can only move along the links 
of a regular lattice defined by the discrete speeds, so that the synchronous particle 
displacements ASi = CiAt never take the fluid particles away from the lattice. For 
the present study, the standard three-dimensional 19-speed lattice is used [1] (see 
Figure [T|) . The right hand side represents the effect of intermolecular solvent-solvent 
collisions, through a relaxation toward local equilibrium, /. , typically a second order 
(low-Mach) expansion in the fluid velocity of a local Maxwellian with speed u: 

(2.12) jf = w iP {l + (3u-c t + ^[uu : {cA - /T 1 ?)]} 

where (3 = mf/ksTf is the inverse fluid temperature (with ks the Boltzmann con- 
stant), Wi a set of weights normalized to unity, and I is the unit tensor in configuration 
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space. The relaxation frequency uj controls the fluid kinematic viscosity v through 
the relation: 

(2.13) v = c 2 s {l/u- At/2) 

where c s is the sound speed in the solvent [7] . Knowledge of the discrete distributions 
fi allows the calculation of the local density p, flow speed pu and momentum-flux 
tensor P, by a direct summation upon all discrete distributions: 

(2.14) P&t)=Y^fi@,t) 

i 

(2.15) pu(x, t) = y] fj(x, t)cj 

i 

(2.16) *P{x,t) =J2fi{x,t)c i c i 

i 

The diagonal component of the momentum-flux tensor gives the fluid pressure, while 
the off-diagonal terms give the shear-stress. Unlike in hydrodynamics, both quantities 
are available locally and at any point in the simulation. 

Thermal fluctuations are included through the source term Fi which reads as 
follows (index notation) 

(2.17) Fi = Wip{F%\ci a c ib - fi-'Sab) + F$ c g iabc } 

where F^ is the fluctuating stress tensor (a 3 x 3 stochastic matrix). Consistency 
with the fluctuation-dissipation theorem at all scales requires the following conditions 

(2.18) (F^(x,t)F%(x' : t')) = 2M^A abcd S(x~x')S(t - if) 

where A a bcd is the fourth-order Kronecker symbol [12]. F^> is related to the fluc- 
tuating heat flux and gi a bc is the corresponding basis in kinetic space, essentially a 
third-order Hermite polynomial (full details are given in [13]). 

The polymer-fluid back reaction is described through the source term Si, which 
represents the momentum input per unit time due to the reaction of the polymer on 
the fluid population ff. 

(2.19) S i (x,t)=Wi(3 iFp+Fp]-^ 

pGD(x) 

where D(x) denotes the mesh cell to which the p th bead belongs. The quantities on 
the left hand side in the above expression have to reside on the lattice nodes, which 
means that the frictional and random forces need to be extrapolated from the particle 
to the grid location. 

The use of a LB solver for the fluid solvent is particularly well suited to this 
problem because of the following reasons: 

i) Free-streaming proceeds along straight trajectories. This is in stark contrast with 
hydrodynamics, in which fluid momentum is transported by its own space-time vary- 
ing velocity field. Besides securing exact conservation of mass and momentum of the 
numerical scheme, this also greatly facilitates the imposition of geometrically complex 
boundary conditions. 
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ii) The pressure field is available locally, with no need of solving any (computationally 
expensive) Poisson problem for the pressure, like in standard hydrodynamics. 

iii) Unlike hydrodynamics, diffusivity is not represented by a second-order differen- 
tial operator, but it emerges instead from the first-order LB relaxation-propagation 
dynamics. The result is that the kinetic scheme can march in time-steps which scale 
linearly, rather than quadratically, with the mesh resolution. This facilitates high- 
resolution down-coupling to atomistic scales. 

iv) Solute-solute interactions preserve their local nature since they are explicitly medi- 
ated by the solvent molecules through direct solvent-solute interactions. As a result, 
the computational cost of hydrodynamic interactions scales only linearly with the 
length of the polymer (no long-range interactions). 

v) Since all interactions are local, the LB scheme is ideally suited to parallel comput- 
ing. 

It is worth mentioning that more advanced Lattice Boltzmann models [141 115) 
could equally well be coupled to the atomistic dynamics. 

2.3. Time exchange. The Molecular-Langevin-Dynamics solver is marched in 
time with a stochastic integrator (due to extra non-conservative and random terms) , 
proceeding at a fraction of the LB time-step, 

dt = At/M 

The time-step ratio M > 1 controls the scale separation between the solvent and 
solute timescales. 

The numerical solution of the stochastic equations is performed by means of a 
modified version of the Langevin Impulse propagation scheme, derived from the as- 
sumption that the systematic forces are constant between consecutive time steps [16 . 
The propagation of the unconstrained dynamics proceeds according to the scheme 
[17| 



, . dt 
r p = r p {t) + —vp{t) 



v* p (t + dt) = e-~< dt !»„(«) + 1 j (F(r p ) + 1Up ) + (> dt / 2 + l) C(dt) 

(2.20) r*(t + dt) = f p + y«p(* + dt) 

where C(dt) is an array of 37V gaussian random variables with zero mean and variance 
^f^-(e 2jdt — 1), and {r} represent temporary positions. The propagator (|2.20p is 
particularly suitable for our purposes since it is second order accurate in time and 
robust, that is, it reduces to the symplectic Verlet algorithm for 7 — > 0. Moreover, at 
variance with the original Langevin Impulse scheme, the modified propagator allows 
for an unambiguous definition of velocities, which are needed to couple the polymer 
to the hydrodynamic field of the surrounding solvent. The particle positions and 
velocities corrected via the SHAKE and RATTLE algorithms read f^(t + dt) — ► f p (t + 
dt) and v£(t + dt) — > v p (t + dt). For consistency, in considering the momentum 
exchange with the solvent the corrected velocities appear in the friction forces. The 
MD cycle is repeated M times, with the hydrodynamic field frozen at time t n = nAt. 

2.4. Spatial exchange. The transfer of spatial information from/to grid to/from 
particle locations is performed at each LB time-stamp t n = nAt. To this purpose, on 
account of its simplicity, a simple nearest grid point (NGP) interpolation scheme is 
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Fig. 2. Transfer of spatial information (a) from grid to particle, and (b) from particle to grid. 
Black spheres denote beads, while in white are the lattice sites. 

used (see FigJ^J). Momentum conservation was checked to hold up to six digits. With 
reference to a time slice t n = nAt, the pseudo-algorithm performing a single LB 
time-step, reads as follows 

1. Interpolation of the velocity: u{x) — » u p 

2. For m = 1,M: 

Advance the molecular state from t to t + dt 

3. Extrapolation of the forces: F p — > F{x) 

4- Advance the Boltzmann populations from t to t + At 
This time-marching can be formally represented by an operator-splitting multi- 
step time procedure for two coupled kinetic equations describing the dynamic evo- 
lution for the fluid and the polymer distribution functions, respectively [TB]. It is 
worth emphasizing that, while LB and MD with Langevin dynamics have been cou- 
pled before, notably in the study of single-polymer dynamics [19], to the best of our 
knowledge, this is the first time such that coupling is put in place for long molecules 
of biological interest. 

2.5. Efficiency considerations. The total cost of the computation scales roughly 

like 

(2.21) t ~ (t LB V + t MD MN)N LB 

where t^s is the CPU time required to update a single LB site per timestep and 
tuD is the CPU time to update a single bead per timestep, V is the volume of the 
computational domain in lattice units and A is the number of polymer beads, with 
M the LB-MD time-step ratio. Finally, Nlb is the number of LB timesteps. In the 
above equation, t^D includes the overhead of LB-MD coupling. Note that tuD is 
largely independent of A because i) the LB-MD coupling is local, ii) the forces are 
short ranged and iii) the SHAKE/RATTLE algorithms are empirically known to scale 
linearly with the number of constraints. 

Regarding the cost of the LB section, this is known to scale linearly with the 
volume occupied by the solvent. For the case where polymer concentration is kept 
constant, the volume needed to accommodate a polymer of A beads should scale 
approximately as A 18 ; however, for translocation studies such as those discussed 
later in this paper, we shall consider a box of given volume, independently on the 
polymer length. 

From the above expression it is clear that M should be chosen as small as possible, 
consistent with the requirement of providing a realistic description of the polymer dy- 
namics. In the present simulation we typically choose M between 5 and 20, depending 
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on the parameters of the simulation, particularly the temperature. This means that 
we are taking the LB representation close to the molecular scale. We will return to 
this important issue in the quantitative discussion of the physical application. 

A tentative estimate of the computational cost proceeds as follows: Assuming 
250 flops/site/LB-step and 2500 nops/bead/MD-step (including the LB-MD coupling 
overhead), and an effective processing speed of 100 Mflop/s, the evolution over 30, 000 
LB steps=150, 000 MD steps of a typical 80 x 40 x 40 grid and 400 beads set-up, would 
take about: 

t = 30, 000 x [250 x (80 x 40 x 4) + 2500 x 5 x 400] /10 8 = (9600 + 1500)sec - 3hrs, 

which is in reasonable agreement with the simulation time observed with the present 
version of the code (7hrs), including the relative MD/LB cost (~ 1 : 4). 

We wish to emphasize that the key feature of the LB-MD approach, namely 
linear scaling of the CPU cost with the number of beads (at constant volume) is 
indeed observed. In fact, the execution times for 50, 100 and 400 beads are 0.433, 
0.489, and 0.882 sec/step, respectively on a 2GHz AMD Optcron processor. By 
excluding hydrodynamics, these numbers become 0.039, 0.075, and 0.318 sec/step. It 
is worth mentioning that thus far, no effort has been directed to code optimization; it 
is quite possible that careful optimization may lower the execution time by an order 
of magnitude. 

2.6. Validation tests. The static and dynamic behavior of the DNA chain 
obtained by our methodology has been compared to the scaling predictions for a single 
chain at infinite dilution. Given the structure factor Sf(k) — Y^i jk e%k ' tyTi ~ r ^), 
standard theory predicts the scaling law Sf(k) — Ng(kR g ), where g(y) is a universal 
function and R g = (1/2N 2 ) j{( r i ~ r j) 2 ) i s the gyration radius. For large k, the 
structure factor is independent of N, and it follows 

S f (k) oc k- 1 ^ 

where experiments, theory and simulations agree on the scaling exponent value /i ~ 
0.584. The static scaling law is not affected by the presence of hydrodynamics. How- 
ever, verification of the scaling law and attainment of the scaling regime for large 
enough chains is a good check for the correctness of our simulation scheme and for 
the subsequent validation of the hydrodynamic behavior. 

The dynamic behavior of the chain is deeply affected by the presence of hy- 
drodynamic interactions. The standard picture of polymer dynamics is based on 
the Rouse (no hydrodynamics) or Zimm (hydrodynamic) description in terms of an 
underlying gaussian chain. In this case, the chain intermediate scattering function 
I s (k,t) = N g~r k \ J2i J -(e* fc 't n w~ ? v( )]) should follow the universal behavior 

I s (k,t)=g(D f R*- 2 tk x ) 

where g(y) is another universal function and Df is the center of mass diffusion con- 
stant. Having introduced the dynamic scaling exponent // via Df ~ , the expo- 
nent x is found to be x = 2 + \x' I \i. According to Zimm theory, fjf = fi and x — 3 [19] 
while, according to Molecular Dynamics simulations, it appears that the actual value 
is somehow lower, i.e. x ~ 2.9 [20] . 

We have considered a chain made of 30 monomers with bead-bead Lennard-Jones 
parameters taken from previous studies of chains simulated via Brownian dynamics 
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Fig. 3. (a) Log-log plot of the structure factor of a polymer in solution made of 30 monomers. 
The straight line is the power law fit in the range 1 < k < 3 with exponent fi = 0.58. (b) Log-log 
parametric plot oft(Ig,k) vs k for 1$ = 0.3 (circles) and Ig = 0.6 (triangles). The lines represent 
the power law fits within the scaling region 1 < k < 3.0 with exponent x = 2.9. 

[2T] (cr = 0.65, e/ksT = 1.0, bond length r = 0.945) in a simulation box of edge 
60. This choice was motivated to verify the range of scaling behavior as compared 
to previous numerical results. We have computed the structure factor, as reported in 
FigGUa), and observed that the scaling regime is clearly visible for 1 < k < 3 with 
an exponent equal /i = 0.58 ± 0.1. In the range of k vectors where the static scaling 
holds, the dynamic scaling has been checked by considering, for given values of the 
computed scattering function, the loci of points t(k,Is) — g~~ 1 {Is)k" x and by fitting 
via a power law curve (see ref. [20] for details). As illustrated in Fig|3jb), the resulting 
scaling exponent is found to be x = 2.9±0.1, in excellent agreement with the expected 
value and similar to previous simulation results on single polymers surrounded by a 
Lattice Boltzmann fluid [21j . Moreover, by applying a heuristic argument [20] . we 
have verified that within the scaling regime, the finite size of the periodic box was not 
biasing the data. 

3. Application: polymer translocation through nanopore. The scheme 
described above is general and applicable to any situation where a long polymer is 
moving in a solvent. This motion is of great interest for a fundamental understanding 
of polymer dynamics in the presence of the solvent. For example, the translocation of a 
polymer through a pore of very small size (of order the separation between monomers) , 
is a process in which the coupling of the molecular motion to the solvent dynamics 
may be of crucial significance. In this section, we will therefore provide a detailed 
discussion of the polymer dynamics in the presence of a solvent for the example of 
translocation through a nanopore but without reference to a specific physical system. 
In the next section we explore the relevance of these results to DNA translocation 
through a nanopore. 

3.1. Initial and Boundary conditions. The polymer is initialized via a stan- 
dard self-avoiding random walk algorithm and further relaxed to equilibrium by stan- 
dard Molecular Dynamics. The solvent is initialized with the equilibrium distribution 
corresponding to a constant density pg and zero macroscopic speed u = 0. 

Boundary conditions for the fluid are periodic at inlet/outlet sections, and zero- 
speed at rigid walls, using the standard bounce-back rule [I]. For the polymer, period- 
icity is again imposed at inlet /outlet, whereas the interaction with rigid walls is han- 
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died by a Lennard- Jones potential with specific wall-polymer parameters a wa ii = 1.5 
and e wa u = 10 -3 in LB units. The connection between slip-flow at the wall and 
intermolecular solid-fluid interactions shall be the objects of future research. 

3.2. Numerical set-up. We consider a three-dimensional box of size N x h x 
N y h x N z h lattice units, with h the spacing between lattice points. We will take 
N x = 2N y , N y = N z ; the separating wall is located in the mid-section of the x 
direction, at x = hN x /2 with N x — 80. At t = the polymer resides entirely in the 
right chamber at x > hN x /2. At the center of the separating wall, a square hole of side 
dhole = 2/i is opened, through which the polymer can translocate from one chamber 
to the other. Translocation is induced by a constant electric field which acts along 
the x direction, and is confined in a rectangular channel of size 2hx h x h along the 
streamline {x direction) and cross-flow (y,z directions). The spatial coarse-graining 
is such that the presence of the solvent as well as electrostatic forces acting due to 
charges on the polymer are neglected altogether as being of secondary importance 
compared to hydrodynamics. 

Here and throughout we work in lattice Boltzmann units, in which length and 
time are measured in units of the lattice spacing h = Ax and time-step At, respec- 
tively. Mass is defined as m = rriLBnisoi- The dimensionless mass tjilb used in the 
simulations is set to unity, which means that mass is measured in units of the solvent 
mass m so i. This choice is not restrictive since the present approach is used to model 
incompressible flows in which density is a parameter which can be rescaled by any 
arbitrary factor. However, it is of some interest to estimate the number of solvent 
molecules represented by a single LB computational molecule, since the inverse of 
this number conveys a measure of the importance of statistical fluctuations at the 
scale of the lattice spacing Ax. Let Sn be this number, which will be defined as 
Psolh — , where plb is the dimensionless density used in the LB simulations. In order 
for the Boltzmann probability distribution to make sense as a statistical observable, 
Nlb » 1- For typical values of p so i = 1 gr/cm 3 , m so ; ~ 20 amu, plb — 1 (which 
correspond to water), and h in the range 1 — 10 2 nm, this yields Sn ~ 10 4 — 10 6 . 
This shows that the neglect of many-body fluctuations inherent to the single-particle 
Boltzmann representation is still justified even at the nanoscopic scale of the lattice 
spacing. 

We will focus here on the fast translocation regime, in which the translocation 
time tx is much smaller than the Zimm time, tz, i.e. the typical relaxation time of 
the polymer towards its native (minimum energy, maximum entropy) configuration. 
Under fast-translocation conditions, the many-body aspects of the polymer dynamics 
cannot be ignored because different beads along the chain do not move indepen- 
dently. As a result, simple one-dimensional Brownian models do not apply [22] . In 
addition to many-body solute-solute interactions, the present approach also takes full 
account of many-body solute-solvent hydrodynamic interactions. The conditions for 
fast-translocation regime can be appraised as follows. The translocation time is esti- 
mated by equating the driving force, F pu u, to the drag force exerted by a solvent with 
dynamic viscosity 77 on a polymer with radius of gyration R, Fd ra g ~ . This 

yields tx ~ Q ZP R . Since the Zimm time is given by tz ~ °, 4 ''^? , the fast-translocation 

"pull f^B -*- 

condition tx <C tz becomes: 

Our reference simulation is performed with F pu u/m — 0.02Aa;/Ai" and ksT/m = 
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10 _ Arc 2 / At 2 , with m the mass of one bead (monomer) of the polymer. The polymer 
length is in the range 20 < N < 400 beads. It can be readily checked that by 
assuming R ~ N - 6 our set of parameters falls safely within the fast translocation 
regime. However, for ksT/rn = 10~ 3 , F pu uR/kBT is of the order of 10 2 — 10 3 which 
is much closer to breaking the above condition. 

The main parameters of the simulation are (in LB units) a — 1.8 and e = 10~ 4 
for the Lennard- Jones potential. The bond length among the beads is set at ro = 1.2. 
According to these values, the Lennard- Jones time-scale, tlj = a/^ ] /2e/m, is of 
the order of ~ lOOAi. Thus, by choosing M — 5 as a time-gap factor, we obtain 
dt ~ txj/500, which is adequate for the resolution of the polymer dynamics. The 
solvent is set at a density p = 1, with a kinematic viscosity v = 0.1Acc 2 /Ai and a 
damping coefficient 7 = 0.1/ At. The flexional rigidity k for the angular potential 
between beads will be 10~ 4 /rad. In order to resolve the structure of the solvent 
accurately on the atomistic scale [53] , we should use a higher resolution of at least 3-4 
orders of magnitude. This means resolving the radial structure of the pore, a task that 
can only be undertaken by resorting to parallel computing. It is nonetheless hoped, 
and verified a posteriori, that this artificial magnification does not affect adversely 
the most significant dynamical and statistical properties of the translocation process, 
by which we mean that eventually, the time-scale of the simulated process may not be 
the same as in the physical process of interest, but the simulated dynamics is related 
to the physical dynamics by a simple rescaling of the time variable. 

3.3. Translocation time. The most immediate quantity of interest in the translo- 
cation process is the dependence of the translocation time on the polymer length. This 
is usually expressed by a scaling law of the form 

T X (N)=t x /dt = N a 

where to is a reference time-scale, formally the translocation time of a single monomer, 
and a a scaling exponent measuring the degree of competition [a > 1) / cooperation 
(a < 1) of the various monomers in the chain. 

We first turn to the derivation of the scaling behavior of the translocation pro- 
cess in the case where hydrodynamic interactions are included. In order to take into 
account the statistical nature of the phenomenon, simulations of a large number of 
translocation events (100 up to 1000) for each polymer length were carried out. The 
ensemble of simulations is generated by different realizations of the initial polymer 
configuration. The duration histograms were constructed by cumulating all events 
for various lengths. Overall, our results are quite similar to the corresponding ex- 
perimental data for DNA translocation through a nanopore [24] . which we discuss in 
more detail in the following section. 

At a next step, our data were shifted and scaled so that the distribution curve 
starts at zero-time and the total probability is equal to unity. The resulting distri- 
butions are on average not gaussians, but skewed towards long translocation times, 
consistent with experiment 24J. Therefore, the translocation time for each length is 
not assigned to the mean time, but to the most probable time, which is the position of 
the maximum in the histogram. In Figf4]the distribution of all the events for polymer 
sizes N=50, 100 and 300 are shown. In this figure, the most probable translocation 
time for each length is denoted by an arrow. From this analysis, a nonlinear relation 
between the most probable translocation time ro and the polymer length is obtained 
that follows closely the theoretically expected scaling tx(N) ~ N a , with a ~ 1.29 
(see Fig [J} . 
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Fig. 4. Probability distributions of the translocation times for various lengths: (a) N =50, 
(b) N=100, and (c) N=300, respectively. Both axes are scaled to produce normalized 
distributions. The arrows show the most probable translocation time for each length. 



100 

number of beads N 



Fig. 5. Translocation time as a function of the number of monomers for the case with hydro- 
dynamics. The straight line represents the power law with exponent 1.29. Time is shown in units 
of the LB timestep At. 



3.4. Dynamics with and without a solvent. A closer inspection into the 
polymer dynamics reveals some interesting features. The molecule shows a blob-like 
conformation on either side of the membrane as it moves through the hole. It may 
either translocate very fast or move from one chamber to the other intermittently, with 
pauses. Both types of events are present with and without a fluid solvent. In addition, 
a careful analysis of all the translocated chains unravels the difference between slower 
and faster translocation within the same fast translocation regime. The nature of the 
variations in time is connected to the random fluctuations of the polymer throughout 
its motion, rather than the temperature or its length. These fluctuations are correlated 
to the entropic forces (gradient of the free energy with respect to a conformational 
order parameter, typically the fraction of translocated beads, see r{t) below) acting 
on both translocated and untranslocated parts of the polymer. In fact, when a solvent 
is present, the interplay between these forces and Fd rag , Fpull determines the motion 
and the shape of the chain and thereby the translocation time. At some point part 
of the chain shapes up in an almost linear conformation increasing in this way the 
entropic force acting on it. This eventually leads to deceleration of the whole chain. 
FigE] shows an illustration of this argument, where a polymer chain, surrounded by a 
solvent is represented at a time where it starts to slow down. In this figure, a polymer 
with the same length but different initial configuration is also shown at the same time. 




Fig. 6. Polymer configuration (N = 400J corresponding to (a) fast and (b) slow translocation 
events. Both snapshots are shown at a tiraestep where the polymer (b) starts to slow down (see 
arrow in Fig\7)). F pu n is applied at the hole region towards a direction indicated by the arrow. 



It is very instructive to monitor the progress in time of the number of translo- 
cated monomers N(i). Note that r(t) = N(t)/N serves as a reaction-coordinate, 
with the translocation time defined by the condition r(tx) = 1. The translocated 
monomers for processes with and without hydrodynamics are shown in Fig|7] For the 
former, events related to the polymers of Fig [6] are shown (curves A\, A3), as well 
as that related to the most probable time (A2). The arrow in this figure indicates 
the timestep corresponding to the snapshots in Fig[5J The translocation for a given 
polymer proceeds along a curve virtually related to its initial configuration and its 
interactions with the fluid. It is clearly visible that there is no general trend. The 
non-hydrodynamic case is in principle different, especially in terms of the time range 
which is larger. This reveals the importance of hydrodynamic coherence. 

Additional insight into the dynamics is obtained by altering the parameter set. 
This has not yet been extensively explored, but it was found that a choice of fc^T = 
10 -5 , Fp U u — 0.01, and e — 0.002 leads to the frequent retraction of the polymer. 
In other words, after having translocated a large fraction of its length, the polymer 
occasionally reverses its motion and anti-translocates away from the hole, never to find 
its way back into it. Moreover, we find that a polymer that retracted in the presence 
of a solvent, manages to fully translocate if the solvent is absent. It is interesting 
to observe that, in principle, no such type of anti-translocating behavior has been 
observed for short polymers. This indicates that hydrodynamics significantly speed-up 
and alter the nature of translocation, especially for long polymers at low temperatures. 
This highly irregular dynamics escapes any scaling or statistical analysis, as well as 
dynamic Monte Carlo simulations |25j . and can only be revealed by self-consistent 
many-body hydro-dynamic simulations. 

4. DNA translocation through a nanopore. The translocation of biopoly- 
mers, such as DNA and RNA plays a major role in many important biological pro- 
cesses, such as viral infection by phages, inter-bacterial DNA transduction, and gene 
therapy [26j . The importance of this process has spawned a number of in vitro exper- 
iments, aimed at exploring the translocation process through micro-fabricated chan- 
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Fig. 7. Progress in time of the number of translocated beads for chains with N = 400 monomers. 
Curves Ai, A3 correspond to slow and fast translocation events (polymers shown in Fig\5^, while 
A2 to an event related to the most probable time. The initial configuration for the polymer in the 
event B is the same as for A2, but in that case no hydrodynamic interactions are included. Time 
is scaled with respect to the value of tq in the case with hydrodynamics. The arrow indicates the 
timestep at which the snapshots in Fig{E[ are shown. 

nels under the effects of an external electric field, or through protein channels across 
cellular membranes [27] • In particular, recent experimental work has focused on the 
possibility of fast DNA-sequencing by reading the base sequence as the polymer passes 
through a nanopore. Some universal features of DNA translocation can be analyzed 
by means of suitably simplified statistical schemes |28j and non-hydrodynamic coarse- 
grained or microscopic models |291 130] . However, a quantitative description of this 
complex phenomenon calls for state-of-the art modeling of the type described above. 
Accordingly, we explore here to what extent the results discussed above for the generic 
situation of polymer translocation apply to the DNA case. 

First, we note that, as already mentioned in the previous section, our results are 
quite similar to the experimental data for DNA translocation through a nanopore 
|24j . Three different interpretations of the current model are physically plausible: 

(a) Following the framework used in recent studies of DNA packing in bacteriophages 
[31] , one monomer in our simulation can be thought of as representing a DNA segment 
of about 8 base-pairs, that is, each bead has a diameter of 2.5 nm, the hydrated 
diameter of B-DNA in physiological conditions. 

(b) It is also physically plausible to assume that a bead represents a portion of DNA 
equivalent to its persistence length of about 50 nm, which translates into mapping 
one bead to ~150 base-pairs. 

(c) Alternatively, as is typically done in simulations of the A-phage DNA in solution 
|32j . one bead can be taken to correspond to ~ 10 3 base-pairs. 

In all three cases h — Ax is equal to the bead size, while the pore, having a width of 
2Ax, will be different from the pores used experimentally, either smaller or larger. In 
addition, the coarse graining model that handles the DNA molecules indicates that 
the MD timescale is stretched over the physical process. A direct comparison between 
our probability distributions for polymer translocation and the experimental results 
sets a different MD timestep for the cases (a), (b), and (c) which is of the order of 3 
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Fig. 8. Translocated (Rt), untranslocated (Rjj) o,nd effective (Re) radii of gyration for the (a) 
non-hydrodynamic and (b) hydrodynamic cases with N = 400. All radii are normalized with respect 
to the initial value Ru(t = 0). Time is also scaled with respect to the total translocation time tx 
for each of the events (a) and (b). The dotted lines denote the regions where Re is nearly constant 
(see text). 



nsec, 100 nsec and 5 /isec, respectively, leading to a LB timestep At = hdt of 15 nsec, 
500 nsec, and 25 /isec. It is difficult at this stage to assign a unique interpretation of 
our model in relation to a physical system. A thorough exploration of the parameter 
space is required before such an assignment can be made convincingly. This is beyond 
the scope of the present work but will be reported in future publications. 

A second encouraging comparison is that the scaling we found for the translo- 
cation time with polymer length (with exponent a = 1.29) is quite close to the 
experimental measurement for DNA translocation [24] (a = 1.27 ±0.03). Beyond 
the apparent consistency between experiment and theory, additional insight can be 
gained by analyzing the polymer dynamics during translocation. 

A hydrodynamic picture of DNA translocation has been presented in ref. [2"4] . 
In this work, the authors assume that the electric held drive is in balance with the 
Stokes drag exerted by the solvent on the blob configuration of the polymer, that is 

Fpuii = D7rpi/ 

T 

where p is the density, v the kinematic viscocity, t the translocation time, and Rjj 
the translocated part of the radius of gyration. In order for this balance to apply at 
all times, it is clear that R 2 must be constant in time, hence it cannot be identified 
neither with the translocated nor with the untraslocated gyration radius of DNA. 
To this end, in Figj8]we represent the time-evolution of the radii of gyration for the 
two sections of DNA, Rtj(t) and -Rt(£), the untranslocated (U) [with x > hN x /2] 
and translocated (T) [with x < hN x /2) parts, respectively. In order to identify a 
time-invariant radius, we define the Rj(t) = QNj~{t), where I = U,T stands for 
the untranslocated and translocated parts and Nj(t) is the corresponding number of 
monomers. The exponent £ ~ 0.6 is the same as previously noted and Q is a constant 
(for long enough polymers). If translocation could be described by the dynamics of a 
single-blob object, characterized by an effective radius of gyration, defined as 



R E (t) = dV<(i), 
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then this quantity should be constant in time. Since the N = Nx(t) + Njj(t) holds 
for all t, the above relations lead to 

(4.1) R E {t) = [i?r /C 0) + i?,y C 0)] C = const 

We focus first on the case without hydrodynamics, FigjS^a). For very small chains 
Ri(t) does not scale as N® 6 and the definition for Re is not valid at the first and 
last ~ 15 — 20% parts of the event, during which the untranslocated and the translo- 
cated parts, respectively, are small. Outside these limits, i?£;(t) as obtained from the 
definition (|4.ip . with the values of Ru(t), Rr(t) directly taken from the simulations, 
is indeed approximately constant. In addition, the values Rr{t =tx) and Ru(t = 0) 
do not coincide, since the former is lower than the latter. This is also the case when 
a solvent is present, FigJS^b). Thus, regardless of the dynamic pathway and the dif- 
ferent conformations the chain may possess during translocation, once the event is 
completed the polymer is more compact than at t = 0. Comparison of the cases with 
and without hydrodynamics reveals that in the latter case the polymer becomes up 
to ~ 7% more confined than when a solvent is added. The untranslocated part of the 
radius of gyration at the end of the process shows an abrupt drop. As a consequence 
the polymer t = tx does not fully recover its initial volume. It it plausible, that 
by allowing the polymer to further advance in time, Rr{t — tx) will become similar 
to Ru(t = to), but this remains to be examined. Nevertheless, in this work we have 
been interested mainly on the chain dynamics related to the first passage times, which 
correspond to the exact period of time needed until all the beads have translocated. 

5. Conclusions. We have presented a multiscale methodology based on the con- 
current coupling of constrained molecular dynamics for the solute biopolymers with a 
lattice Boltzmann treatment of solvent dynamics. Owing to the dual field-particle na- 
ture of the Lattice Boltzmann technique, this coupling proceeds seamlessy in time and 
only requires standard interpolation/extrapolation for information-transfer in phys- 
ical space. This multiscale methodology has been applied to the case of polymer 
translocation through a nanopore, with special emphasis on the role of hydrodynamic 
coherence on the dynamic and statistical properties of the translocation process. It is 
found that hydrodynamic interactions play a major role in accelerating the translo- 
cation process, especially for long molecules at low temperature. 

An attempt to connect these results to the process of DNA translocation through 
a nanopore revealed certain similarities with experiment, especially in the scaling 
law of the translocation time with polymer length. The presence of hydrodynamic 
interactions lead to a decrease in the translocation times, compared to the cases 
without a fluid solvent. Inspection of the variation of the translocated beads and the 
radii of gyration with time reveals interesting aspects of the DNA dynamics during 
translocation. 

Future directions for the simulations include the detailed study of the effects of 
temperature, finite-length and geometrical details of the nanopore geometry, as well as 
electrostatic interactions of the DNA molecule with the surrounding fluid. To this end, 
resort to parallel computing is mandatory, and we expect the favourable properties of 
LB towards parallel implementations to greatly facilitate the task. Work along these 
lines is currently in progress. 
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